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Robust Neighboring Optimal Guidance 
for the Advanced Launch System 

David G. Hull* 

In recent years, optimization has become an engineering tool through 
the availability of numerous successful nonlinear programming codes. 
Optimal control problems are converted into parameter optimization 

(nonlinear programming) problems by assuming the control to be 
piecewise linear, making the unknowns the nodes or junction points of the 
linear control segments. Once the optimal piecewise linear control 
(suboptimal) control is known, a guidance law for operating near the 

suboptimal path is the neighboring optimal piecewise linear control 
(neighboring suboptimal control). Research conducted under this grant has 
been directed toward the investigation of neighboring suboptimal control 
as a guidance scheme for an advanced launch system. The list of 
references is a list of papers presented at technical meetings; these papers 
are included at the end of the report. 

The first step is to obtain the optimal piecewise linear control for the 
advanced launch system, upon which the neighboring optimal piecewise 

linear control is based. These results have been obtained by using a 

nonlinear programming code based on recursive quadratic programming 
with numerical partial derivatives and are reported in Ref. 1. In an effort 
to improve the results obtained by nonlinear programming, the suboptimal 
control problem is solved by the shooting method which requires 
analytical derivatives (Ref. 2). By guessing the control history, the 
shooting method is completely desensitized to the guesses of the initial 
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Lagrange multipliers. Ref. 2 shows that the optimal piecewise linear control 
can be a good approximation of the actual optimal control. 

Once the suboptimal control has been obtained, the next step is to 
develop the neighboring suboptimal control for guidance about the 
suboptimal path. Since this is a completely new research area, a simpler 
launch model, the lunar launch problem, has been used to determine the 
feasibility of this guidance strategy. 

First results are reported in Ref. 3. Given a perturbation in the state 
from the state corresponding to the suboptimal control, the control 
parameter perturbations are obtained by minimizing the increase in the 
performance index (second variation) subject to the constraint that the 
final conditions are satisfied. This process leads to a set of gains which 
multiply the state perturbations to get the control parameter 
perturbations. At this point, time is the running variable, and while the 
results are satisfactory, they give some indication that the method has not 
yet been applied properly. 

Refs. 4 and 5 are further steps to clarify the fundamental issues of 
this guidance approach. In Ref. 5, the problem is reduced to a “fixed final 
time” problem by using the horizontal component of velocity as the 
variable of integration. These results indicate that neighboring suboptimal 
control is now being formulated and applied correctly. 

Continuing work is associated with using the time as the running 
variable since this is the way most guidance systems operate. 

Finally, a list of master’s degrees awarded during this research effort 
is given in Table 1. 
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Abstract 

The maximum-final-mass trajectory of a proposed 
configuration of the Advanced Launch System is pre- 
sented. A model for the two-stage rocket is given; the 
optimal control problem is formulated as a parame- 
ter optimization problem; and the optimal trajectory 
is computed using a nonlinear programming code called 
VF02AD. Numerical results are presented for the con- 
trols (angle of attack and velocity roll angle) and the 
states. After the initial rotation, the angle of attack goes 
to a positive value to keep the trajectory as high as pos- 
sible, returns to near zero to pass through the transonic 
regime and satisfy the dynamic pressure constraint, re- 
turns to a psotive value to keep the trajectory high and 
to take advantage of minimum drag at positive angle of 
attack due to aerodynamic shading of the booster, and 
then rolls off to negative values to satisfy the constraints. 
Because the engines cannot be throttled, the maximum 
dynamic pressure occurs at a single point; there is no 
maximum dynamic pressure subarc. 

To test approximations for obtaining analytical solu- 
tions for guidance, two additional optimal trajectories 
are computed: one using untrimmed aerodynamics and 
one using no atmospheric effects except for the dynamic 
pressure constraint. It is concluded that untrimmed 
aerodynamics has a negligible effect on the optimal tra- 
jectory and that approximate optimal controls should 
be able to be obtained by treating atmospheric effects 
as perturbations. 

List of Symbols 

Cl lift coefficient 
Co drag coefficient 
Cm pitching moment coefficient 
D aerodynamic drag (lb) 
g local gravitational acceleration (ft/sec 2 ) 
h altitude (ft) 
i orbit inclination (rad) 

I tp vacuum specific impulse (sec) 

J performance index 
l aerodynamic reference length (ft) 
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L aerodynamic lift (lb) 
m vehicle mass (slugs) 

Ma aerodynamic pitching moment (ft lb) 

M Mach number 
P penalty function 
p atmospheric pressure (lb/ft 2 ) 
q dynamic pressure (lb/ft 2 ) 

Sb aerodynamic reference area (ft 2 ) 

T thrust (lb) 

T va t vacuum thrust (lb) 
t $ staging time (sec) 

V velocity (ft/sec) 

a angle of attack (rad) 

7 flight path angle (rad) 

6 thrust gimbal angle (rad) 

0 pitch angle (rad) 

A longitude (rad) 
p velocity roll angle (rad) 
p atmospheric density (slug/ft 3 ) 
r latitude (rad) 
r normalized time 

w rotational velocity of earth (rad/sec) 
heading angle (rad) 

Subscripts 

b body axes 
eg center of gravity 
e exit 
/ final 
i inertial 
o initial 
8 sea-level 
w wind axes 

I. Introduction 

A program is under way to develope an unmanned, 
all-weather, launch system for placing medium to large 
payloads (^ 120,0001b) into low-earth orbit. A prospec- 
tive design for this Advanced Launch System (ALS) is 
shown in Fig. 1 to be composed of a core vehicle and 
a booster. Both the booster and the core are ignited 
at launch, and staging occurs when all the booster pro- 
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pellant is consumed. Payload mass can be increased by 
adding another booster. 

Part of the design process is to iterate the vehicle de- 
sign and trajectory design until a reasonable combina- 
tion is achieved. This paper is concerned solely with the 
optimal trajectory design of the proposed configuration. 
The objective is to find the trajectory which maximizes 
the final mass (since the engines burn throughout the 
trajectory, this is also a minimum final time problem). 
Any remaining propellant can be considered for conver- 
sion to payload or a decrease in launch weight. The 
physical model is that of flight over a rotating, spher- 
ical earth with an exponential atmosphere. Launched 
vertically from the surface of the earth, the payload is 
to be placed into perigee of an 80nm by 150nm transfer 
orbit. Because of structural considerations, there is a 
limit on the amount of dynamic pressure the vehicle can 
withstand. 

This study has had several goals: (a) to determine the 
maximum-final-mass trajectory of the proposed ALS, 

(b) to generate initial Lagrange multipliers for a shooting 
code to investigate neighboring extremal guidance, and 

(c) to determine if atmospheric effects (pressure thrust 
and aerodynamics) can be considered as a perturbation 
to vacuum thrust and gravity for guidance law devel- 
opment. While only (a) and (c) are reported here, (b) 
requires the use of an exponential atmosphere. Hence, 
the dynamic pressure limit based on a standard atmo- 
sphere has been lowered to have the same effect in an 
exponential atmosphere. 

In Section 2, a model is presented for the proposed 
ALS configuration. Then, the optimal control problem 
is formulated in Section 3 and converted into a param- 
eter optimization problem in Section 4. This is done 
for relative ease in obtaining an optimal trajectory. Nu- 
merical results are presented in Section 5 in the form 
of optimal controls, states, and dynamic pressure. Also 
contained in Section 5 are two additional optimal trajec- 
tories based on untrimmed aerodynamics and neglected 
atmospheric effects. Finally, conclusions are presented 
in Section 6 . 


II. Physical Model 

In this section, a physical model for the Advanced 
Launch System (ALS) is defined. It includes the equa- 
tions of motion for flight over a rotating, spherical earth 
with an exponential atmosphere and the mass, propul- 
sion, and aerodynamic properties of the vehicle. 

Equations of Motion 

Since sideslip causes drag, the vehicle is assumed to fly 
at zero sideslip angle, so that only the angle of attack 
gives the orientation of the vehicle relative to the free 
stream. The direction of the lift vector is then controlled 


through the bank angle or, more specifically, through the 
velocity roll angle. 

A three-degree-of- freedom model for vehicle motion 
can be obtained from a six- degree-of- freedom model by 
one of two aerodynamic approximations: untrimmed 
aerodynamics or trimmed aerodynamics. For a rocket, 
untrimmed aerodynamics is equivalent to setting the 
thrust gimbal angle to zero and ignoring the aerody- 
namic pitching moment. On the other hand, with 
trimmed aerodynamics, it is assumed that the pitch rate 
is zero (pitching moment equals zero) so that the gim- 
bal angle can be determined as a function of the angle 
of attack. 

In view of the above comments, the three-degree-of- 
freedom equations of motion relative to the earth are 
gjven by (Ref. 3) 
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In these equations, A is the longitude, r is the latitude, h 
is the altitude above mean sea level, V is the velocity, 7 
is the flight path angle, t/> is the heading angle, m is the 
mass, r = r, + h is the distance from the center of the 
earth to the vehicle center of gravity, w is the angular 
velocity of the eath, D is the drag, L is the lift, T is the 
thrust, I, p is the specific impulse, 6 is the gimbal angle 
of the thrust vector, a is the angle of attack, and /i is the 
velocity roll angle. With regard to signs, a positive roll 
angle generates a negative heading toward the south. 

For trimmed aerodynamics, the pitching moment, 
which is the sum of the aerodynamic pitching moment 
and the thrust pitching moment, is set equal to zero, and 
the resulting expression solved for the thrust gimbal an- 
gle. With reference to Fig. 1 and by assuming that 6 is 
small, this process leads to 
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Atmosphere 



Figure 1: Force and Moment Nomenclature 


where Ma is the aerodynamic pitching moment ant l t 
is the distance from the center of gravity to the exit 
plane of the engines. Because 6 is dependent on the 
aerodynamic pitching moment and the moment is de- 
pendent on the pitching moment coefficient, it results 
that 6 is linear in a with the coefficients varying with 
time. Aerodynamics is discussed in further detail later 
in this section. 

Eqs. (1) have two singularities: V s 0 in the 7 and 
the ij> equations and 7 = y in the equation. To re- 
move the V singularity and to clear the launch tower, 
the vehicle is flown vertically for 3 sec with the angle of 
attack and the bank angle being chosen so that 7 = 0 
and = 0. To remove the 7 singularity, the vehicle is 
pitched over at constant heading = 0 ) for 1.0 sec at 
a constant negative pitch rate 6 whose optimal value is 
determined. Since 6 = 7 + a, the angle of attack during 
pitch-over is given by 

a = | - 7 + <>(<- 3) . (3) 

Finally, the bank angle is chosen to make = 0. With 
a flat earth model, p = 0 . 

Earth 

The earth is taken to be a rotating, spherical body 
whose surface is described by the mean sear level radius 
r, and whose gravitational acceleration varies with alti- 
tude according to the inverse-square law 

> = ( 4 ) 

where g,r j represents the earth’s gravitational parame- 
ter. Sea-level gravitational acceleration g,, r, , and the 
rotational velocity of earth u are known constants given 

“ ft 

r, = 2.09256725E+7 ft , g, = 32.174 — ? 

* sec^ 

w = 7.2921 158E-5 — . (5) 


The atmosphere is represented by the exponential 
functions 


— = exp(-^) , ?- = exp(-£) 

P, A, p, "2 


( 6 ) 


where the scale-height constants are given by 

Ai = 23,800 ft , A 2 = 23,200 ft (7) 

and the searlevel values of the density and pressure are 

( 8 ) 

it' 

Finally, the speed of sound is given by 


p, = .002377^ , p. =2, 116.24 ^ . 


D 

a =^ ~p 


( 9 ) 


where 7 = 1.4 is the ratio of specific heats of air. 


Mass Characteristics 

The ALS configuration consists of a core vehicle as 
depicted in Fig. 1. The take-off mass of the ALS con- 
sists of the inert vehicle mass, the propellant mass, pay- 
load mass, payload margin mass, and the payload fairing 
mass (Table 1). 


Table 1: Mass Characteristics 


Vehicle 

Vehicle Component 

Take-off Mass 
(slugs) 

Core 

Inert Mass 
Propellant 
Payload 

Payload Margin 
Payload Fairing 

5,474.29 

45,974.38 

3,729.71 

372.97 

1,215.89 


Total Core 

56,767.26 

Booster 

Inert Mass 
Propellant 

6,740.85 

45,066.82 


‘'Total Booster 

51,807.67 

Core + 
Booster 

Total Take-off Mass 

108,574.93 


The center of gravity is located relative to a coordinate 
system whose origin is at the tip of the core vehicle, 
whose x axis is down the longitudinal axis, and whose 
y axis is toward the booster. For the first stage, the 
vehicle center of gravity is assumed to have coordinates 

x cy = 165.45 ft , y cg = 10.36 - .0388* ft (10) 
so that l t is constant and has the value 

l t = /-x c , = 110.81 ft (11) 


sec 
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where / = 276.26 ft is the length of the core vehicle. 
Actually, x cg varies slightly but this variation has been 
neglected. For the second stage, untrimmed aerodynam- 
ics is used so the eg position is not needed. 

Propulsion 

The ALS is powered by ten liquid hydrogen/liquid 
oxygen low cost rocket engines (LCE): seven power the 
booster and three power the core. All engines are ig- 
nited at launch; staging occurs when the booster fuel is 
depleted; and the core engines bum until insertion. 

Propulsion characteristics of interest are the thrust 
T, vacuum thrust T vac , and the specific impulse 7, p (see 
Eqs. 1). If the exit pressure is conservatively approxi- 
mated as p t = 0, the thrust of a single engine is modeled 
as 

T'=T W „'-M/ (12) 

where the prime denotes one engine, p is the atmospheric 
pressure at the altitude of the rocket, and A / is the exit 
area. Date relevant to one LCE are as follows: 

T va c= 580, 110.0 lb 
A.' = 40.381 ft 2 (13) 

I t p = 430.0 sec . 


C m = C m (A/,a) (16) 


where M denotes the Mach number and the bar indi- 
cates that the moment is about a fixed point (launch 
eg). About the actual center of gravity, the moment is 


given by 


Cm — Cd 


Vcg - 10.36 

l 


(17) 


since x tg is assumed not to change. 

While the aerodynamic data could have been used in 
tabular form with linear interpolation to read the tables, 
the approach taken is to assume polynomials in a with 
Mach- number-dependent coefficients. For the first stage, 
the coefficients are written as 


Cp = Cd 0 (M) + Cz) fta (M)a 2 + Cd^{M)q 
Cl = Ci a (M)ot (18) 

^m=^m 0 (A/) + ^m.(M)a 


where the Mach-number-dependent terms have been ob- 
tained from cubic-spline curve fits of the tabular data. 
After staging, the flow regime is hypersonic and the aero- 
dynamic force coefficients are modeled as 

Cd = Cdo + Cd + Cd ^ q7 

Cl = Cl* ® d - Cl** o? (19) 


For the complete vehicle, 

T = kT* , I$p sss I t p , T va c = kT V ae (14) 

where Jfc = 10 before staging and k = 3 after staging. 
Specific impulse is like specific propellant consumption 
(weight flow rate of propellant per pound of thrust); 
hence, it has the same value regardless of the number 
of engines operating. 


where the coefficients of a are constants. Also, pitching 
moments are assumed to be negligible after staging, that 
is, untrimmed aerodynamics are used (£ = 0). 

A peculiarity of the aerodynamics of the combined 
vehicle at supersonic and hypersonic speeds is that the 
drag coefficient has a minimum at a positive angle of 
attack. This is caused by the aerodynamic shading of 
the booster by the flow field of the core. 


Aerodynamics 

The drag, lift, and pitching moment are related to 
their respective coefficients by the standard equations 

D = qSbCD > L = qSbCL , A I a = qSblC m (15) 

where q = \pV 2 is the dynamic pressure, S\, = 
1413.71 ft 2 is the cross-sectional area of the combined 
vehicle (booster 4- core), and / is the length of the core. 
While the aerodynamic coefficients are needed at and 
about the center of gravity (eg), the aerodynamic data 
has been provided at and about the launch eg. Although 
the drag and lift transfer directly, the moment changes 
with eg position. Therefore, the aerodynamic data at 
the eg must be related to the launch eg. 

The aerodynamic data are preliminary estimates as- 
sociated with the development of the six-degreee-of- 
freedom simulation presented in Ref. 4. These data are 
provided in tabular form (Tables 2 through 6) consistent 
with the functional relations 

Cp = C|?(Af,a) } Cl = Cl( Af,a) , 


III. The Optimal Control Problem 

Formally the optimal control problem considered here 
is to find the control history u(t) which minimizes a per- 
formance index of the form 

J = $(xj) (20) 

subject to the differential constraints 

i = /(x,u), (21) 

the prescribed initial conditions 

to — to, , x 0 — x 0 , » (22) 

the prescribed final conditions 

$(x ; ) = 0 , (23) 

and a state-variable inequality constraint 

S(x) < 0 . (24) 
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Each of these quantities is discussed below. 


State Variables and Control Varables 

The state variables are x T = [A r h V y rp m] while 
the control variables are u T = [a /i]. 

Performance Index 


It is desired to maximize the final mass. Hence, the 
performance index is taken to be 


$ = - 


rrirej 


(25) 


where the minus sign is included because the perfor- 
mance index is actually minimized and where m re y is 
the sum of the payload mass, the payload margin mass, 
and the payload fairing mass. A performance index of 
<$ = — 1.0 means that the reference mass is inserted into 
orbit with no extra fuel. 


Differential Constraints 


The differential constraints are the equations of mo- 
tion (Eqs. 1) completely expressed in terms of the state 
variables and the control variables. 

Prescribed Initial Conditions 

For the trajectory design problem, the initial condi- 
tions are taken to be 

t 0 = Osec , A 0 = -80.54 deg ,r 0 = 28.5 deg 

h a = Oft ,V 0 = 0 — , 7 <> = 90.0 deg (26) 
sec 

= 0.0 deg ,m 0 = 108, 574.93 slugs 

During the vertical rise segment, the heading angle is 
undefined, so the initial condition on \p is actually the 
heading angle during the pitch-over segment. 


Prescribed Final Conditions 


The Advanced Launch System is being designed to 
place a nominal payload at perigee of an 80nm by 150nm 
transfer orbit of 28.5 deg inclination. As a consequence, 
the equality constraint residuals are 

= h, -486,080 ft , = V h -25,776.9 — 

* sec 

4*3 = yj , = cosij - cos(28.5deg) (27) 

where the inertial velocity and the inclination are related 
to the relative states as follows (Ref. 5 and 6) 

y. + 2V rucosycosxpcosr + (rwcosr) 2 ] ^ (28) 


cosi = 


cos r ( V cosy cos + rwcosr) 

[V 2 cos 2 y -f 2Vru> cosy cosxp cost -f ( tujcost ) 2 ] i 


St ate- Variable Inequality Constraint 

Based on structural considerations, the ALS must 
not exceed a maximum dynamic pressure of q = 
650 lb/ft 2 . Therefore, the state- variable inequality con- 
straint residue 5 is 

S = ^pV 2 - 650 lb/ft 2 . (29) 

Actually, in a standard atmosphere, the limit is q ma x = 
850 lb/ft 2 . The value of 650 lb/ft 2 is chosen because the 
value of p is approximately 20% smaller in the exponen- 
tial atmosphere than the standard atmosphere around 
the maximum dynamic pressure portion of the trajec- 
tory. 


IV, The Suboptimal Control Problem 


The optimal control problem is converted to a param- 
eter optimization problem (suboptimal control problem) 
as follow: (a) the time is normalized by introducing the 
transformation r = (b) the control ti(f) is replaced 

by a set of nodal points which is linearly interpolated, 
and (c) the state-variable inequality constraint is con- 
verted to a point constraint by using a penalty function. 

Because of the time transformation, the boundary val- 
ues of r are given by 

A 3 4 

To ~ 0 , Tp — — - , T\ — , 

tj if 


153.54 

r, = — , Tf s 1 

l S 


(30) 


where t p = 3 sec is the time at the beginning of pitch- 
over and t\ ss 4 sec is the time when three-dimensional 
flight begins. Staging occurs when all of the booster 
propellant is consumed; hence, t $ = 153.54 sec. 

Figure 2 shows the arrangement of nodal points in 
each stage. Nine nodes are used for the control during 
the first stage, and five for the control during the second 
stage. Even though the duration of the first stage is 
shorter than that of the second, there is more activity in 
o during the first stage, making more nodes desireable. 
The nodes are equally spaced in each stage so that the 
node times are 


Ti 


n + 



- 1) , i = 1 


9 


n = t. + 1 — -(»- io) , i = io- 14 . (3i) 

4 


5 



u (rad) 



Figure 2: Example Control History 



(38) 


( 39 ) 


Derivatives required by the nonlinear programming al- 
gorithm are computed by central differences. 


Note that there are two control nodes at the stage time. 
This has been done in order to find the true suboptimal 
control. 

The dynamic pressure constraint is converted to a pa- 
rameter inequality constraint by introducing the penalty 
function 

P = — f min 2 [(l — ),0]<ft>0 (32) 

Jt 9 

which accumulates value when q > q m & *• The constraint 
becomes 

P f > 0 . (33) 

To compute Pj , the penalty function is differentiated to 

form a 

P = — min 2 [(l ),0] (34) 

Qmas 

where 

Po = o . (35) 

In all, the nonlinear programming problem involves 30 
parameters, that is, the parameter vector is given by 

X = [(9 , ai , . . . , «i4 , i • • • i . «/] (36) 

where 6 is the pitch rate during the pitch-over, a* , m 
are the angle of attack and the bank angle nodes, and 
tj is the final time. 

If values of the parameters (36) are known, the differ- 
ential equations (1) and (34) can be integrated through 
the mission to determine the states and P at the fi- 
nal time. Then, the performance index (25), the orbital 
insertion equality constraint residuals (27), and the dy- 
namic pressure inequality constraint (29) can be com- 
puted. It follows that the performance index and the 
constraints are functions of the parameters (36) such 
that the nonlinear programming problem can be ex- 
pressed as follows: 

Find the set of parameters X which minimizes the per- 
formance index 

j _ _ m /,W (37) 

m re j 

subject to the equality constraints 

h J. 


V. Numerical Results 

The optimal trajectory has been computed using a 
nonlinear programming code known as VF02AD which is 
based on quadratic programming. Optimal control his- 
tories are presented in Fig. 3, while the resulting states 
are shown in Figs. 4 through 7. The magnitude of the 
performance index is 103.94% where 100% = 171,120 
lb. This means that an additional 6,742 lb of payload 
can be placed in orbit with this vehicle by using the op- 
timal trajectory. The vehicle is inserted into orbit at 
tj = 363.8 sec and the optimal value of the pitch rate 
during the 1.0 sec pitch-over is -.02005 rad/sec. 

Shown in Fig. 8 is the dynamic pressure. It is seen 
that the maximum dynamic pressure occurs at a single 
point and not along a q m ax subarc. This is due to the 
no-throttling design of the vehicle and the fact that the 
aerodynamic forces needed to fly along q = 9mar cannot 
be achieved. Optimal trajectories with lower values of 
q mmB have been calculated, and the results are the same. 

It is difficult to completely determine the meanings of 
the optimal control histories because performance-index 
minimization and constraint satisfaction are going on 
all through the trajectory. For angle of attack, it is seen 
from Fig. 3 that the vehicle initially goes to positive a to 
achieve altitude and decrease q. Then, the dip in a from 
t = 40 to 60 sec allows the vehicle to pass through the 
transonic regime efficiently (Mach 1 occurs at t « 50 sec) 
and to satisfy the dynamic pressure inequality constraint 
(flmax occurs at t « 70 sec). Next, the vehicle returns to 
positive o to get low drag and to decrease the magnitude 
of 7. Staging occurs around Mach 8 and the roll off in a 
from positive to negative values during the second stage 
helps pull the trajectory down to meet the final condi- 
tions. For the velocity roll angle, the nonzero values at 
the beginning of the trajectory seem to be caused by the 
rotational effects of earth where the vehicle wants to fly 
at constant latitude throughout most of the first stage. 
Changes in /1 near the end of the trajectory help cause 
constraint satisfaction, particularly in the orbit inclina- 
tion. 

Additional optimal trajectories have been computed 
with the intent of determining what kinds of approxi- 
mations can be made in order to obtain approximate 
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analytical solutions for guidance purposes. First, the ef- 
fect of using untrimmed aerodynamics (6 = 0) rather 
than trimmed aerodynamics is shown in Fig. 9 and 10 to 
change only slightly the optimal controls and to cause a 
relative change in the performance index of 0.2% (376.5 
lb). Hence, untrimmed aerodynamics is a reasonable 
approximation. Second, the question of whether or not 
atmospheric effects can be considered a perturbation is 
considered. This means that the pressure term in the 
thrust and the aerodynamics are neglected; however, the 
dynamic pressure constraint is maintained because it is a 
structural constraint. The optimal controls for this case 
are shown in Fig. 1 1 and 12 and lead to a relative increase 
in the performance index of 16% (27,379 lb). Trajectory 
profiles for the atmosphere and no-atmosphere cases are 
shown in Fig. 13. The optimal control which results 
from the no-atmosphere case is reasonably close to that 
of the atmosphere case and has the same general trend. 
This seems to indicate that atmospheric effects can be 
treated as a perturbation. 

VI. Discussion and Conclusions 

The maximum-final-mass trajectory has been com- 
puted for a two-stage rocket representing the Advanced 
Launch System and operating over & rotating, spherical 
earth with an exponential atmosphere. The problem is 
converted into a parameter optimization problem by re- 
placing the control histories by node points and using 
straight-line interpolation to form functions. Then, a 
nonlinear programming code known as VF02AD is used 
to perform the optimization. Optimal trajectories have 
been calculated for three cases: (a) trimmed aerody- 
namics, (b) untrimmed aerodynamics, and (c) no atmo- 
sphere. With the assumption of trimmed aerodynamics, 
the aerodynamic model is as accurate as possible for a 
three-degree-of-freedom analysis. The optimal trajec- 
tory is characterized by positive angles of attack over 
most of the path with a prominant decrease during pas- 
sage through maximum dynamic pressure. The maxi- 
mum dynamic pressure occurs at a single point rather 
than over a sub arc because the engines cannot be throt- 
tled. 

To obtain analytical solutions for guidance purposes, 
approximations must be introduced. The effect of re- 
placing trimmed aerodynamics by untrimmed aerody- 
namics has been examined, and it is concluded that 
untrimmed aerodynamics gives good results. 

Next, the effect of neglecting atmospheric effects 
(pressure thrust and aerodynamics) has been investi- 
gated. With the exception of the transonic and max- 
imum dynamic pressure portion of the trajectory, it is 
clear that atmospheric effects can be considered as per- 
turbations to the trajectory generated by vacuum thrust 
and gravity. During the passage through the transonic 
and maximum dynamic pressure part of the trajectory, 


there is a difference of 14 deg between the atmosphere 
and no-atmosphere solutions. Since this region consti- 
tutes less than fifteen percent of the whole trajectory, 
treating atmospheric effects as perturbations could yield 
satisfactory results. 
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Figure 3: TYimmed Aerodynamic Control Histories 
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Altitude (ft) 
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Figure 5: Flight Path and Heading Angle vs. Time; 
Trimmed Flight 



Time (sec) 

Figure 6: Altitude and Velocity vs. Time; Trimmed 
Flight 
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Figure 7: Latitude and Longitude vs. Time; Trimmed 
Flight 



Figure 8: Dynamic Pressure vs. Time; Trimmed Flight 



Figure 9: Angle of Attack vs. Time 
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Figure 10: Velocity Roll Angle vs. Time 
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Figure 13: Trajectory Profiles; Atmosphere and 
No-atmosphere 



Time (sec) 


Figure 11: Angle of Attack vs. Time 



Table 2. Lift Coefficient (core + booster) 


Subsonic Data 


Angle of Attadc (deg) 


M 

±0.0 

±2.0 

±4.0 

±6.0 

±8.0 

0.0 

o.6 


0.1775“ 

0.2663 

0.355 

0.2 

0.0 

0.08876 

0.1775 

0.2663 

0.355 

0.4 

0.0 

0.08876 

0.1775 

0.2663 

0355 

0.6 

0.0 

0.08876 

0.1775 

0.2663 

0355 

0.8 

0.0 

0.08876 

0.1775 

0.2663 

0.355 

1.0 

0.0 

0.08720 

0.1744 

0.2616 

0.3488 


Supersocuc/Hypersonic Data 


± 10.0 

0.4438 

0.4438 

0.4438 

0.4438 

0.4438 

0.4360 


M 

0.0 

Angle of Attack (deg) 
2.0 4.0 6.0 

8.0 

10.0 

1.2 

0.0 

0.0862 

0.1724 

0.2586 

03448 

0.431 

1.5 

0.0 

0.086 

0.171 

0.260 

0351 

0.431 

2.0 

0.0 

0.090 

0.175 

0.262 

0.354 

0.435 

2.5 

0.0 

0.098 

0.181 

0.268 

0370 

0.460 

3.0 

0.0 

0.100 

0.192 

0.278 

0.385 

0.490 

3.5 

0.0 

0.102 

0.200 

0.290 

0.401 

0.510 

4.0 

0.0 

0.104 

0.202 

0.291 

0.405 

0.510 

5.0 

0.0 

0.104 

0.206 

0.298 

0.410 

0.509 

6.0 

0.0 

0.103 

0.203 

0.300 

0,408 

0.508 

7.0 

0.0 

0.100 

0.195 

0.298 

0.400 

0.502 

8.0 

0.0 

0.095 

0.185 

0.290 

0.395 




Angle of Attack (deg) 



M 

-2.0 

-4.0 

-6.0 

-8.0 

-10.0 



1.2 

1.5 
2.0 

2.5 

3.0 

3.5 

4.0 

5.0 

6.0 

7.0 

8.0 


-0.084 

-0.086 

-0.090 

-0.098 

- 0.100 

- 0.120 

- 0.120 

- 0.120 

-0.125 

-0.115 

- 0.110 


-0.170 

-0.171 

-0.175 

-0.181 

-0.192 

- 0.200 

-0.215 

-0.225 

-0.225 

- 0.222 

-0.218 


-0.260 

-0.260 

-0.262 

-0.268 

-0.278 

-0.290 

-0.310 

-0.327 

-0.334 

-0.332 

-0.325 


-0.350 

-0.351 

-0.354 

-0.370 

-0385 

-0.401 

-0.420 

-0.442 

-0.451 

-0.452 

-0.450 


-0.431 

-0.431 

-0.435 

-0.460 

-0.490 

-0.510 

-0.520 

-0.542 

-0.567 

-0.580 

-0.565 


Figure 12: Velocity Roll Angle vs. Time 
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Table 3. Drag Coefficient (core + booster) 


M 


- 2.0 


Angle of Attack (deg) 
- 4.0 - 6.0 


- 8.0 


Subsonic Data 


M 


±0.0 


Angle of Attack (deg) 
±2.0 ±4.0 ±6.0 


±8.0 


0.0 

0.2 

0.4 

0.6 

0.8 

1.0 


0.1870 

0.1872 

0.2062 

0.2599 

0.3480 

0.7800 


0.1904 

0.1906 

0.2096 

0.2633 

0.3514 

0.7834 


0.2024 

0.2026 

0.2216 

0.2753 

0.3634 

0.7954 


0.2254 

0.2256 

0.2446 

0.2983 

0.3864 

0.8184 


0.262 

0.2622 

0.2812 

0.3349 

0.4230 

0.8550 


±10.0 

0.314 

0.3142 

0.3340 

0.3877 

0.4758 

0.9078 


Supersonic/ Hypertonic Data 


M 


0.0 


Angle 

2.0 


of Attack (deg) 
4.0 6.0 


8.0 


10.0 


0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

1.2 

1.5 
2.0 

2.5 

3.0 

3.5 

4.0 

5.0 

6.0 

7.0 

8.0 


- 0.0201 

- 0.020 

-0.019 

-0.018 

-0.016 

-0.004 

0.003 

0.009 

0.009 

0.007 

0.005 

0.004 

0.004 

0.005 

0.008 

0.008 

0.008 


-0.044 

-0.044 

-0.043 

-0.042 

-0.040 

-0.027 

-0.029 

-0.019 

-0.0155 

-0.016 

-0.018 

-0.018 

-0.019 

-0.018 

-0.017 

-0.017 

-0.017 


- 0.067 

- 0.067 

- 0.067 

-0.066 

-0.064 

-0.051 

-0.058 

-0.048 

-0.045 

-0.043 

-0.041 

-0.040 

-0.040 

-0.038 

-0.028 

-0.028 

-0.028 


- 0.091 

-0,091 

- 0.091 

-0.089 

-0.087 

-0.075 

-0.089 

-0.077 

-0.071 

-0.067 

-0.063 

-0.062 

-0.062 

-0.058 

-0.058 

-0.058 

-0.058 


- 10.0 

- 0.115 


- 0.115 

- 0.115 

- 0.113 

- 0.111 

-0.098 

-0.119 

-0.106 

-0.097 

-0.092 

-0.089 

-0.086 

-0.085 

-0.082 

-0.078 

-0.076 

-0.075 


1.2 0.800 

0.805 

0.815 

0.838 

0.875 

0.928 






1.5 0.740 

0.703 

0.645 

0.640 

0.635 

0.635 


Table 5. 

Lift Coefficient (core 

vehicle) 

2.0 0.672 

0.656 

0.555 

0.525 

0.525 

0.525 






2.5 0.648 

0.628 

0.512 

0.468 

0.465 

0.455 



Angie of Attack (deg) 



3.0 0.637 

0.608 

0.486 

0.448 

0.431 

0.418 

M 

±0.0 

±2.0 ±4.0 ±6.0 

±8.0 

±10.0 

3.5 0.630 

0.596 

0.470 

0.425 

0.406 

0.392 

8.0 

0.2062 

0.2089 0.2206 0.2417 

0.2733 

0.3160 

4.0 0.628 

0.587 

0.460 

0.410 

0.385 

0.368 

10.0 

0.2180 

0.2201 0.2313 0.2523 

0.2835 

0.3262 

5.0 0.620 

0.572 

0.448 

0.392 

0.355 

0.352 

12.0 

0.2353 

0.2374 0.2486 0.2703 

0.3024 

0.3459 

6.0 0,617 

0.570 

0.446 

0.382 

0.348 

0.348 






7.0 0.615 

0.567 

0.445 

0.378 

0.340 

0.340 

Table 6. 

Draff Coefficient (core vehicle) 

8.0 0.615 

0.565 

0.445 

0.372 

0.340 

0.338 














Angle of Attack (deg) 



Angle of Attack (deg) 



M 

±0.0 

±2.0 ±4.0 ±6.0 

±8.0 

±10.0 

M -2.0 

-4.0 

-6.0 

-8.0 

-10.0 


8.0 

0.2062 

0.2089 0.2206 0.2417 

0.2733 

0.3160 

1.2 0.803 

0.815 

0.838 

0.875 

0.928 


10.0 

0.2180 

0.2201 0.2313 0.2523 

0.2835 

0.3262 

1.5 0.745 

0.750 

0.773 

0.800 

0.871 


12.0 

0.2353 

0.2374 0.2486 0.2703 

0.3024 

0.3459 

2.0 0.690 

0.708 

0.731 

0.768 

0.822 







2.5 0.665 

0,680 

0.706 

0.745 

0.790 







3.0 0.648 

0,651 

0.688 

0.730 

0.771 







3.5 0.640 

0.650 

0.675 

0.716 

0.757 







4.0 0.631 

0.641 

0.665 

0.706 

0.745 







5.0 0.625 

0.635 

0.651 

0.692 

0.731 







6.0 0.610 

0.625 

0.646 

0.686 

0.727 







7.0 0.610 

0.620 

0.640 

0.685 

0.730 







8.0 0.610 

0.620 

0.640 

0.684 

0.725 








Table 4. Pitching Moment Coefficient (core + 
booster) 


Angle of Attack (deg) 


M 

0.0 

2.0 

4.0 

6.0 

0.0 

0.0035 

0.0271 

0.0508 

0.0744 

0.2 

0.0035 

0.0271 

0.0508 

0.0744 

0.4 

0.0040 

0.0276 

0.0513 

0.0745 

0.6 

0.0052 

0.0288 

0.0538 

0.0757 

0.8 

0.0072 

0.0308 

0.0558 

0.0777 

1.0 

0.020 

0.046 

0.072 

0.098 

1.2 

0.033 

0.062 

0.093 

0.123 

1.5 

0.038 

0.066 

0.095 

0.124 

2.0 

0.033 

0.059 

0.085 

0.111 

2.5 

0.030 

0.056 

0.077 

0.097 

3.0 

0.029 

0.052 

0.071 

0.087 

3.5 

0.028 

0.049 

0.066 

0.080 

4.0 

0.027 

0.047 

0.061 

0.076 

5.0 

0.026 

0.045 

0.057 

0.068 

6.0 

0.026 

0.042 

0.054 

0.068 

7.0 

0.0255 

0.042 

0.053 

0.068 

8.0 

0.0255 

0.042 

0.052 

0.068 


8.0 10.0 


0.0981 

0.1217 

0.0981 

0.1217 

0.0986 

0.1222 

0.0998 

0.1234 

0.1018 

0.1254 

0.124 

0.150 

0.153 

0.183 

0.153 

0.182 

0.134 

0.162 

0.120 

0.135 

0.103 

0.116 

0.094 

0.099 

0.0895 

0.099 

0.085 

0.098 

0.082 

0.096 

0.082 

0.096 

0.082 

0.097 
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A SHOOTING APPROACH TO SUBOPTIMAL CONTROL 

David G. Hull' an<l Jyh-Jong Shorn* 

Department of Aerospace F.nginooring aivl Engineering Mechanics 
The University of Texas at Austin 


Abstract 

The shooting method is used to solve the suboptimal control prob* 
lcm where the control history is assumed to be piecewise linear. 
Suboptimal solutions can be obtained without difficulty and can 
by increasing the number of nodes lead to accurate approximate 
controls and good starting multipliers for the regular shooting 
method. Optimal planar launch trajectories are presented for the 
Advanced Launch System. 


Formally, this fixed-final-time suboptimal control problem is 

to find the parameters 1/ which minimize the perfor- 

mance index 

J m #*/, t/) + £ U T > x, u„ . . . ,m«, */)dr (4) 

subject to the dynamics 

/(r, W 


1. Introduction 

The original motivation for using the shooting method to solve the 
suboptimal control problem (piecewise linear control) has been 
to calculate an accurate suboptimal control and ultimately to 
find the corresponding neighboring extremal feedback control rule. 
Since aerospace minima are usually quite flat, an approximate op- 
timal control can deliver most of the optimal performance. Then, 
the ability to compute the suboptimal control and the neighboring 
extremal without difficulty would be useful. 

In this paper, the shooting method is developed for the subop- 
timal control problem and used to optimize the Advanced Launch 
System trajectory. The usual sensitivity of the solution process 
to the initial guesses disappears completely, and solutions are ob- 
tained without difficulty. Of course, only an approximate optimal 
control is achieved, but if it is not good enought, its accuracy can 
be improved by increasing the number of control nodes. 

2. Suboptimal Control Problem 

The standard optimal control problem is to find the control u(f) 
which minimizes the Scalar performance index 

J m <>(*/,</) + r *.«)<& G) 

subject to the system dynamics 


(7) 


the prescribed boundary conditions 

TfesO, zq * So,, r/« 1, “ °- 

In these equations, the prime denote, a derivative with reepect to 

T ’“ 4 Kw, u m tj) - 

/(r,x,U| u.,1/) - <//(«/’•.*♦“) 

where 

«(T)-«a + “ W — (r-n), n<rsn*. (8) 

and the node time* r* ere fixed. 

By the usual arguments of the calculus of variations, the equa- 
tions defining the suboptimal solution are given by 


✓ -/ 
v.-BJ 

« 0 

p fttjdr « -G,, G 


/f-I + AU 


(9) 


and 


r« ■ 0, So * *o.. T J ’ 


i, • o, V* 


( 10 ) 


x*s/(<,*,u), (2) 

and the prescribed boundary condition* 

to - 0, *0 ■ *«. . ¥»(*/.</) ■ 0 • W 

The dimensions of z, u, and r!> are n X 1, r X 1, and p X 1, respec- 
tively. This problem is made into a suboptimal control problem 
by normalizing the final time through the transformation r ■ t/t f 
and by restricting the class of functions to which the optimal con- 
trol can belong. Here, the restricted class is that of piecewise 
linear functions The end points U| of the straight line 
segments are called nodes. 

l MJ. Thomp*on Regent* ProfcMor 
^Graduate Rcaeaxch AwiataDt 


I. Shooting Matbod 

ro put the suboptimal control problem in a form suitable for ap- 
plying the shooting method, new states v*(r) and u>(r)are in- 
troduced to eliminate the integrals in Eq. (9). The optimality 
conditions become 


tnd 


**■/. 
vi - 

w' w Hu 


k m l,...,m 


T 0 l 

T / 


= 0, 
1 1, 


*8 : 

Vk, 


* * 0 ., 

0, 

* 0, 


v*. 

x r 

w, 


. 0, 

*-Gu 


Wg ■ 0 


(ID 


( 12 ) 


1 



5 Conclusions 


If a new stale vector :(r) is d. iim d as 

:* „ [r T \ T v, .... H-D 

ami a parameter vector is introduced as 

n T = (ui ... u m </] , 

the differential equations (II) can be rewritten in the form 

z' = F(t,z,o) ( 15 ) 

Of the initial states, there are 1 + n + m + 1 conditions; only 
A« is unknown. At the final time, there are in Eqs. (12) 1 + P + 
n + m + 1 final conditions. Of these, p equation* are solved for 
the p Lagrange multipliers i/ which are in turn eliminated from 
the remaining conditions to form 

h(z,,a)* 0 ( 16 ) 

whose dimension is (n + m + 1) x 1. . 

The derivation of the equations for the shooting method is 
straightforward and leads to the following algorithm: 

1. Guess Ao and a 

2. Integrate from r 0 = 0 to tj = 1 

2 * s F z o known 

<t>j = F,*i * 2 , = to 1 0 0) T (17) 

<I>' = F,# + F. *o ** 0 

3. Calculate ||A||. 

4. Calculate S\ 0 and 6a by solving 

(MV *•/*/! [ “ ~ ah ° 8) 

and using a norm reduction scheme to determine a. 

5. Check for convergence (||h|| < e)- If not, go to 2. 

The advantage of this method is that there is absolutely no 
influence of Ao on *. On the other hand, the sensitivity of the 
shooting method to Ao is replaced by having to accept an approx- 
imate solution. However, by using a reasonable number of nodes, 
it should be possible to obtain Ao’s for which the exact shooting 
method can be converged. 


4. Optimal Planar Trajectory for the ALS 

The Advanced Launch System is a two-stage rocket consisting of 
a core with a side-mounted booster. Staging occurs at the fixed 
time of burnout of the booster. Ref. 1 contains a description of 
the physical model. 

In the optimization problem, the performance index is the final 
mass; the state equations are the equations of motion for flight in 
a great circle plane over a nonrotating spherical earth 
control is the angle of attack; the initial conditions are all 
and final conditions are imposed on altitude, velocity, and flight 
path angle* 

Converged results are presented in Table l and Fig. 1 for three 
first and second stage node arrangements. Starting multipliers for 
the 3-2 case are given in Table 1, and the control nodes are taken 
to be Q = -.5, 10., 6., 4., -4. deg and l, = 300. sec. Convergence 
required 17 iterations and 241 sec of CPU time on a CDC Cyber 
computer. Also presented are the converged values obtained from 
the standard shooting method. Note that the optimal results are 
approached as the number of nodes is increased. 

Other than having to derive the multiplier equations, no dilh- 
cutties have been encountered during these calculations. 


The shooting approach to suboptima) control is an effective way 
to obtain approximate optimal trajectoiies and to obtain starting 
Lagrange multipliers for the regular shooting method. 
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Table 1: Converged Results 



node pattern _ j 

optimum 

3-2 suets 

3-2 

5-5 


p.l 


0.8522 

6.4545“ 

0.8541 

0.8544 

^(sec^ 

300.00 

371.72 

371.69 

371.64 i 

3 ( i .C>3 


1.0 

-8.662E-6 

-8.925fc-6 

-8.914E-IT 

•8 .939E-G 

*Xo 

1.0 

-4.148E-T 

-4.150E-4 

*4129£-T 

-4.125E-4 

A* 

1.0 

1. C6lL.il" 

8.702E-3 

3.925E-i 

2.584E-3 

ya — 
^mO 

1.0 

•2.004E-5 

-2.017E-5 

-2.0166-5 

*2.0 1 3E-5 



Figure 1; Angle of Attack Histories 
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ABSTRACT 

The neighboring extremal feedback control law is de- 
veloped for systems with a piecewise linear control for 
the case where the optimal control is obtained by non- 
linear programming techniques. To develop the control 
perturbation for a given deviation from the nominal 
path, the second variation is minimized subject to the 
constraint that the final conditions be satisfied. This 
process leads to a feedback relationship between the 
control perturbation and the measured deviation from 
the nominal state. A simple example, the lunar launch 
problem, is used to demonstrate the validity of the guid- 
ance law. For model errors on the order of 5%, the 
results indicate that 5% errors occur in the final condi- 
tions. 

INTRODUCTION 

In order to develop the neighboring optimal guidance 
law for a dynamical system, it is first necessary to ob- 
tain the optimal control, and this can be a formidable 
task. Currently, most trajectory optimization is ac- 
complished by restricting the class of control functions 
to some subclass, say piecewise linear functions (sub- 
optimal control). Then, the control variables are pa- 
rameters (nodes of piecewise linear function), and the 
suboptimal control is found by applying nonlinear pro- 
gramming methods. Hence, the subject of this paper 
is the development of the neighboring suboptimal feed- 
back control law, assuming that the suboptimal control 
law is available. 

Given the suboptimal control and a perturbation in 
the state at some time, the neighboring suboptimal con- 
trol is found by minimizing the increase is the perfor- 
mance index subject to the constraint that the final 
conditions must be satisfied. Since the first variation 
vanishes, minimizing the increase in the performance 
index is equivalent to minimizing the second variation. 
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The constraint of satisfying the final conditions is ob- 
tained through the use of transition matrices to the final 
point. The above process leads to an analytical expres- 
sion for the gains of the neighboring suboptimal feed- 
back control law. Because of the simplicity of the con- 
trol law, the suboptimal control rule can be applied to 
the vehicles rather than sample and hold. This should 
allow the sample time to be increased, if errors do not 
grow too rapidly. 

To test this guidance rule, it is applied to a simple 
trajectory problem with various levels of modeling er- 
rors. The results indicate that this guidance approach 
has merit. 

RTTBOPTIMAL CONTRO L PROBLEM 

The optimal control problem [1] being considered 
here is to find the control history u(t) which minimizes 
the performance index 

J = <(>(tj,Xf) (1) 

subject to the state differential equations 

i = /(<,*,«), ( 2 ) 

the prescribed initial conditions 

to = to,. *0 = * 0 . , (3) 

and the prescribed final conditions 

V’(t/,*/) = o. ( 4 ) 

Here, this problem is converted into a suboptimal con- 
trol problem [2] by assuming that the controls are piece- 
wise linear, meaning that the unknowns become the 
junction points (nodes) of the linear control segments 
and the final time. 

If a denotes the unknown parameter vector, that is, 
a T = [tj, tin . «i2, ••• . «ai. « 22 , ••• ], the suboptimal 
control problem is stated as follows: 

Find the set of parameters a which minimizes the 
performance index 

J = F(a) (•'>) 
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subject to the equality constraints 

C(a ) = 0 . (6) 

The differential constraints are an integral part of defin- 
ing the functions F and C and are written as 

dx / \ 

- = 9 {r,x,a) ( 7 ) 

To =0, X 0 = X0.» T J - 1 

where r = i/tf and x 0 . are the specified values of the 
initial states. 

It is assumed that this problem is solved numerically 
by using a nonlinear programming code, and the next 
step is to find the neighboring suboptimal feedback con- 
trol law. 

NEIGHBORING SUBOPTIMAL CONTROL 

The solution of the suboptimal control problem gives 
nominal control and state histories to be followed by the 
vehicle. However, because of modelling errors, the ve- 
hicle when using the nominal control deviates from the 
nominal state. Hence, it is desired to find the neighbor- 
ing suboptimal control perturbation which enables the 
vehicle to operate in the neighborhood of the nominal 
trajectory. The general philosophy is to find the con- 
trol perturbation which minimizes the increase in the 
performance index while satisfying the prescribed final 
conditions. 

Since the first variation vanishes along the subopti- 
mal path, the increase in the performance index is the 
second variation 


A J = ±6a T G aa 6o (8) 


where G — $ + u r 'i is the augmented performance in- 
dex and v is a constant Lagrange multiplier. Once the 
suboptimal control has been obtained numerically, the 
second derivative matrix G a „ can be computed numer- 
ically. The next step is to find the constraints on 6a 
which guarantee satisfaction of the final conditions (4). 

The variation of the state equation (7) leads to the 
differential equation 


— 6x = g x 6x - f g a 6a 
dr 


( 9 ) 


which must be solved subject to the boundary condi- 
tions 


TO — To, , 
TJ = 1, 


Sxq — &Zq 9 
\l> x ,bxj -f IpTjbtf 


= 0 . 


Next, the solution of Eq. (9) is assumed to have the 
transition matrix form 


bx — $bxj + 

(11) 

where 

<P/ = /, V/ = o 

(12) 

to guarantee that Sxj = 6xj . Then, substituting Eq. 
(11) into Eq. (9) and equating like coefficients leads to 
the following differential equations: 

4 = g x $ 

¥ = + 9a 

(13) 

which must be solved subject to the boundary condi- 
tions (12). Once $ and have been obtained, Eq. (11) 
can be used. 

To satisfy the final condition (10), Eq. (11) is rewrit- 
ten as 

6xj = $~ l 6z — (14) 

Then, assuming 'J’, / = 0, Eq. (10) leads to 


if> Xj $~ 1 6x — tl> ltf Q~ l 'i6a = 0 . 

(15) 

Applied to ro, this equation becomes 


V>, y *o ~ 1>s,*o l6 *o = 0 

(16) 


and is the constraint on the control node perturbation 
6 a imposed by the final condition. 

The last step is to minimize A J as given by Eq. (8) 
with respect to 6a subject to the constraint (16). Stan- 
dard parameter optimization methods lead to 

6a = K 0 6x 0 (1*0 


where the gain Kq is given by 

K 0 = Gj a 1 *j4>o T V>?’/ 

(V-*, *o 1 *<>G - a 1 *0*0 1 * I , )' : 1 ^*i *0 1 • 


(18) 


If the sampling is performed continuously, the param- 
eter perturbation becomes 


6a = K 6x 


(19) 


where 


K 


Gji' • 

( tfir , 1 1/£, )- 1 


( 20 ) 


These gains can be computed at several values of r 
and stored in the onboard computer for interpolation 


( 10 ) 
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purposes. 



Two difficulties occur in the use of Eq. (19) as a 
guidance law. First, * goes to zero as r approaches 
unity so that the computation of the gains becomes 
indefinite (zero over zero). This has been handled in 
the following application by computing the gains at r = 
.950 and r = .975 and extrapolating them to r = 1. 
The second problem is determining the value of r on 
the perturbed path since the perturbed final time is 
unknown. This has been accomplished iteratively by 
guessing 6tj, computing t = f/(f/ + ^/)> computing 
6a and, hence, 6tj , and repeating the computation until 
the computed 6tj nearly equals the guessed 6tf. 

EX AMPLE - LUNAR LAUNCH PROBLEM 

The lunar launch problem has been selected as a sim- 
ple example to illustrate the application of this guidance 
law. The optimal control problem is to find the thrust 
inclination history 0(t) which minimizes the time to in- 
sertion 

J — tj (21) 

subject to the differential constraints 


1 a 
l(^) 



% Deviation from Optimal j| 


Integrate 

Control 

19.760 

-5.0 

y 

msm 

1.15 



u 

5.96 

6.35 



V 

0.88 

0.81 J 

1 21.840 

+5.0 

m 

0.99 

0.96 




im 

4.45 

iigg| 


V 

0.34 

1.00 


Table 1: Results for 5% Modeling Error in Thrust 
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% Change 

State 

% Deviation from Optimal 

{■&*) 

in g 


Sample- 

Hold 

Integrate 

Control 

5.054 

- 5.0 

y 

0.07 

0.08 



u 

0.11 

0.26 



V 

0.25 

0.01 

5.586 

+5.0 

y 

0.10 

0.09 



u 

0.43 

0.41 



V 

0.17 

0.07 


Table 2: Results for 5% Modeling Error in Gravity 


x — ti 

y = v 

u = a cos 9 
v — a sin 9 — g , 

the prescribed initial conditions 


t 0 = x 0 = Po = «o = vo = 0 , 


( 22 ) 


(23) 


and the prescribed final conditions 

y 0 = 50, 000 ft, u 0 = 5,444 ft/sec, v 0 = 0 ft/sec. 

(24) 

The quantities a and g are the constant thrust acceler- 
ation and lunar acceleration of gravity whose nominal 
values are a — 20.8 ft/sec^ and g — 5.32 ft/sec . 

Using five nodes for the suboptimal control calcula- 
tion leads to 

tj = 272.7 sec, 

f?i = 26.09 deg, 

$2 = 20.68 deg (251 

$ 3 = 15.34 deg, ' 1 

04 ss 9.061 deg, 

$3 = 3.113 deg 

To test the guidance law, a 5% error is introduced 
in a which drives the vehicle away from the nominal. 
Gains have been computed stored at each .025 in r. 
Two implementations have been performed: one is to 


use sample and hold and the other is to use the actual 
linear control. Results are shown for a 4 sec sample 
time in Table 1. Note that a 5% error in a leads to 
roughly a 5% error in the insertion conditions. 

That the linear control does not do uniformly better 
t h An sample snd hold is disappointing. It is felt that 
the sample time could be increased substantially for the 
linear control relative to sample and hold and still yield 
good results. At any rate these are preliminary results 
and further study is warranted. 

Similar results have been developed for a 5% error in 
g and are shown in Table 2. Qualitatively, these results 
are similar to those in Table 1. 


pjfiriTRSTON AND CONCLUSIONS 

The neighboring extremal feedback control law has 
been developed for systems with a piecewise-linear con- 
trol whose nominal control and trajectory have been 
computed using nonlinear programming techniques. 
Given a perturbation in the state, the neighboring ex- 
tremal control perturbation is obtained by minimizing 
the increase in the performance index relative to the 
nominal value subject to the constraint that the final 
conditions be satisfied. Numerical results for the lunar 
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launch problem with mismatches in the thrust accel- 
eration and gravity acceleration show that 5% model 
errors lead to 5% final condition errors. Further study 
of this guidance law seems warranted. 
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Abstract 

The neighboring extremal feedback control law is 
developed for systems with a piecewise linear control 
for the case where the optimal control is obtained by 
nonlinear programming techniques. To develop the 
control perturbation for a given deviation from the 
nominal path, the second variation is minimised sub- 
ject to the constraint that the final conditions be sat- 
isfied. This process leads to a feedback relationship 
between the control perturbation and the measured 
deviation from the nominal state. 

Introduction 

In order to develop the neighboring optimal guid- 
ance law for a dynamical system, it is first neces- 
sary to obtain the optimal control. Currently, most 
trajectory optimisation (see Ref. 1 for example) is 
accomp lish ed by restricting the class of control func- 
tions to some subclass, say piecewise linear functions 
(suboptimal control). Then, the control variables are 
parameters (nodes of piecewise linear function), and 
the suboptimal control is found by applying nonlin- 
ear programming methods. Hence, the subject of this 
paper is the development of the neighboring subop- 
timal feedback control law, assuming that the sub- 
optimal control law is available. 

Suboptimal Control Problem 

The optimal control problem being considered here 
is to find the control history u(r) which minimises 
the performance index 

J = (1) 

subject to the state differential equations 

= /(r, *,«,</) , (2) 

the prescribed initial conditions 

T 0 = To,, *0 = *0., (®) 
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and the prescribed final conditions 

tj - 1 , rp(x f ,tj) = 0. ( 4 ) 

Here, the time has been normalised by the final time, 
that is, r = t/tf where tj is an unknown parame- 
ter. This optimal control problem is converted into a 
suboptimal control problem (parameter optimisation 
problem) by assuming that controls are piecewise lin- 
ear, meaning that the unknowns become the nodes 
of the linear control segments and the final time. 

If a denotes the unknown parameter vector, that 
is, a T = [f/, 1 * 11 , * 13 , ••• , *31 > > ••• ], differ- 

ential equations (2) and its boundary conditions can 
be rewritten as 

^ = §(r,x,a), to = t 0 ., *o = *o., Tj . 
iT ( 6 ) 
Given «, these equations can be integrated to obtain 
xj = */(«) so that d = d[*/(«),*/l = *X a ) 
d = y»[*/(a),</] = C(a) Then, the suboptimnl con- 
trol problem is to find the parameter vector a which 
minimises the performance index J = F(a) subject 
to the constraint C(«) = 0. 

To solve the suboptimal control problem ana- 
lytically, the augmented performance index J ' — 
F(a) + v T C(a) = G(a, */) is formed. The first vari- 
ation conditions are G, = 0 and C — 0 which de- 
termine a and v. The second variation becomes 
S 7 J’ = 6a T G a .6a > 0 where C.6a = 0. 6a can be 
divided into dependent and independent parts, and 
the second variation condition becomes the positive 
definiteness of a matrix. 

At this point, it is assumed that the suboptimal 
control problem is solved by using a nonlinear pro- 
gramming code (see Ref. 1 , for example), and the 
next step is to find the neighboring suboptimal con- 
trol. 

Neighboring Suboptimal Control 

The solution of the suboptimal control problem 
gives nominal control and state histories to be fol- 
lowed by the vehicle. However, because of modelling 
errors, the vehicle when using the nominal control 



deviates from the nominal state. Hence, it is desired 
to find the neighboring suboptimal control pertur- 
bation which enables the vehicle to operate in the 
neighborhood of the nominal trajectory. The gen- 
eral philosophy is to find the control perturbation 
which minimizes the increase in the performance in- 
dex while satisfying the prescribed final conditions. 

Since the first variation vanishes along the subop- 
timal path, the increase in the performance index is 
the second variation 

AJ = \ 6a T G aa 6a (6) 


and is the constraint on the control node perturba- 
tion 6a imposed by the final condition. 

The last step is to minimize AJ as given by Eq. 
(6) with respect to 6a subject to the constraint (14). 
Standard parameter optimization methods lead to 

6a — K$6xo (15) 

where the gain Kq is given by 

K. = G;. l *3Wr,- 


subject to C a 6a = 0 which is imposed below. Once 
the suboptimal control has been obtained, the second 
derivative matrix G« a can be computed numerically. 
The next step is to find the constraints on 6a which 
guarantee satisfaction of the final conditions (4). 

The variation of the state equation (5) leads to the 
differential equation 

~ 6z = g x 6x + g a 6a (7) 

dr 

which must be solved subject to the boundary con- 
ditions 

Tb = ro, , 6x o = 6zq b 

Tj = 1, 1>*,6xj + i>t s 6tf = 0 . 

Next, the solution of Eq. (7) is assumed to have the 
transition matrix form 

6z = *6zf + 96a (9) 


where 

*/ = I, t/ = 0 (10) 

to guarantee that 6zj = 6zj. Then, substituting Eq. 
(9) into Eq. (7) and equating like coefficients leads 
to the following differential equations: 


*' = 9x * 

*' = + g e 


( 11 ) 


which must be solved subject to the boundary con- 
ditions (10). Once + and * have been obtained, Eq. 
(9) can be used. 

To satisfy the final condition (8), Eq. (9) is rewrit- 
ten as 

= ( 12 ) 

Then, for the case where ip tj = 0, Eq. (8) leads to 
xP t/ *- l 6z-rp tJ *-'*6a = 0. (13) 

Applied to tq, this equation becomes 

4>r J *t''*o6a-1>r,*o l 6z o = 0 (14) 


Application 

In Ref. 2, neighboring suboptimal control has been 
applied in the same manner as neighboring optimal 
control, that is, sampling is assumed to occur contin- 
uously so that ro = r. However, in optimal control, 
any part of an optimal trajectory to the final con- 
straint manifold is an optimal trajectory, but this is 
not the case in suboptimal control. In fact, there 
may not even be enough nodes between the sample 
point and the final constraint manifold to satisfy the 
boundary conditions. 

Two alternate approaches are being considered. 
First, additional nodes are placed near the final 
constraint manifold to make neighboring suboptimal 
control valid near the end of the trajectory. Second, 
the suboptimal control is computed from each node 
to the final constraint manifold, and the gains (16) 
are computed at each node. These gains are linearly 
interpolated far the operation of the vehicle. Unfor- 
tunately, no results for either case are available at 
the time of this writing. 
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The neighboring optimal feedback control law is developed for sys- 
tems with a piecewise linear control for the case where the optimal 
control is obtained by nonlinear programming techniques. To de- 
velop the control perturbation for a given deviation from the nomi- 
nal path, the second variation is minimized subject to the constraint 
that the final conditions be satisfied (neighboring suboptimal con- 
trol). This process leads to a feedback relationship between the 
control perturbation and the measured deviation from the nomi- 
nal state. Neighboring suboptimal control is applied to the lunar 
launch problem. Two approaches, single optimization and multiple 
optimization, for calculating the gains are used, and the gains are 
tested in a guidance simulation with a mismatch in the acceleration 
of gravity. Both approaches give acceptable results, but multiple 
optimization keeps the perturbed path closer to the nominal path* 

INTRODUCTION 

In order to develop the neighboring optimal guidance law for a dynamical system, 
it is first necessary to obtain the optimal control. Currently, most trajectory opti- 
mization (see Ref. 1, for example) is accomplished by restricting the class of control 
functions to some subclass, say piecewise linear functions (suboptimal control). Then, 
the control parameters are the nodes of a piecewise linear function, and the subop- 
timal control is found by applying nonlinear programming methods. The subject of 
this paper is neighboring optimal control for systems with piecewise linear controls, 
or neighboring suboptimal control, and its application to vehicle guidance. 

In Refs. 2 and 3, the neighboring suboptimal control problem is formulated as a 
free final time problem and applied to the lunar launch problem. This formulation 
requires an iteration at each sample point to find the normalized time. In this paper, 
neighboring suboptimal control is formulated as a fixed final time problem and applied 
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to the lunar launch problem. While this problem is a minimum time problem, it can 
be converted to a “fixed final time” problem by using the horizontal component of 
velocity, whose final value is fixed, as the variable of integration. 

Two approaches for computing the control gains are presented. In the single 
optimization approach, the nominal suboptimal control is viewed as a sequence of 
reduced-node suboptimal controls to the final constraint manifold. Hence, the quality 
of the suboptimal control diminishes along the flight path. In the multiple optimiza- 
tion approach, a new full-node suboptimal control is computed from each node of the 
nominal suboptimal trajectory to the final constraint manifold. Hence, the quality of 
the suboptimal control along the flight path is maintained. 

After the suboptimal control problem and the neighboring suboptimal control 
problem are summarized, the lunar launch problem is defined. Then, the single 
optimization and multiple optimization approaches are used to compute the gains 
which are, in turn, tested in a simulation with a mismatch in the acceleration of 
gravity. Finally, some conclusions are reached about the use of these two approaches. 

SUBOPTIMAL CONTROL PROBLEM 

The fixed final time optimal control problem being considered here is to find the 
control history u(r) which minimizes the performance index 

J = <f>{xj) (!) 

subject to the state differential equations 


Tr = 


the prescribed initial conditions 


to = t 0 ,, x 0 = xo„ v°; 

and the prescribed final conditions 

Tf = 1, Tp{x f ) = 0. ( 4 ) 

Here, the time has been normalized by the final time, that is, r = t/tj. This op- 
timal control problem is converted into a suboptimal control problem (parameter 
optimization problem) by assuming that controls are piecewise linear, meaning that 
the unknowns become the nodes of the linear control segments. 

If a denotes the unknown parameter vector which for one control is written as 
a T = [tii, u r ], the differential equation (2) and its boundary conditions can be 

rewritten as . 

dx t \ 

— = g(T,x,a) 

To = T 0j , X 0 = X 0j , Tf = 1. 
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( 5 ) 

( 6 ) 



Given a, these equations can be integrated to obtain i/ = xj(a) so that J = 
<^[a; / (a)] = F(a) and ip[xj(a)} = C(a) Then, the suboptimal control problem is to 
find the parameter vector a which minimizes the performance index J = F(a ) sub- 
ject to the constraint C(a ) = 0. 

To solve the suboptimal control problem analytically, the augmented performance 
index J' = F(a) + v T C(a ) = G{a, //) is formed. The first variation conditions are 
G a = 0 and C = 0 which determine a and v. The second variation becomes 8 2 J' = 
Sa T G a a^ a > 0 where C a 6a = 0. 8a can be divided into dependent and independent 
parts; the dependent parts can be eliminated; and the second variation condition 
becomes the positive definiteness of a matrix. 

At this point, it is assumed that the suboptimal control problem is solved by using 
a nonlinear programming code (see Ref. 1, for example), and the next step is to find 
the neighboring suboptimal control. 

NEIGHBORING SUBOPTIMAL CONTROL 

The solution of the suboptimal control problem gives nominal control and state 
histories to be followed by the vehicle. However, because of modeling errors, the 
vehicle deviates from the nominal state. Hence, it is desired to find the neighboring 
suboptimal control perturbation which enables the vehicle to operate in the neigh- 
borhood of the nominal trajectory. The general philosophy is to find the control 
perturbation which minimizes the increase in the performance index while satisfying 
the prescribed final conditions. 

Since the first variation vanishes along the suboptimal path, the increase in the 
performance index is the second variation 

A J = \8a T G aa 8a (7) 

where the second derivative matrix G aa can be computed numerically. The elements of 
8a are not independent but are constrained by the need to satisfy the final conditions 

6xl> = rJ> Xj 6x f = 0. (8) 

The variation of the state equation (5) leads to the differential equation 

-7- 8x = g x 8x + g a 8a (9) 

dr 

which must be solved subject to the boundary conditions 

To = T 0l , 8x 0 = 8x 0, (10) 

tj = 1 , = 0 - 

Next, the solution of Eq. (9) is assumed to have the transition matrix form 

8x = Q>8x j + ty8a (11) 
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where 


$; = /, 'iff = 0 . (12) 

to guarantee that Sxj = Sxj. Then, substituting Eq. (11) into Eq. (9) and equating 
like coefficients leads to the differential equations 

= g x if+g a ( 13 ) 

= 9x$ 

which must be solved subject to the boundary conditions (12). Once $ and ^ have 
been obtained, Eq. (11) can be used. 

To satisfy the final condition (10), Eq. (11) is evaluated at To and rewritten as 

Sxj = $ 0^0 - ( 14 ) 

Then, Eq. (10) leads to 

^x ) ^o 1 ^oSa-ip x/ i>o 1 Sx o = 0 ( 15 ) 

which is the constraint on the control node perturbation, 8a , imposed by the final 
condition. 

The last step is to minimize A J as given by Eq. (7) with respect to 8a subject to 
the constraint (15). Standard parameter optimization methods lead to 

8a = K 0 8xq (16) 

where the gain Kq is given by 

Ko = (17) 

The computation of the gains can be checked by observing that K 0 = da^/dzo 
and using numerical differentiation. Given a suboptimal control and state history, 
a perturbation in the state is introduced at some node, and the suboptimal control 
from that perturbed state to the final constraint manifold is computed. The gains 
are computed as K 0 (i,j) = Aa(i)/Ax 0 (j) where Aa is the change in the suboptimal 

control caused by the change in the state. 

The application of neighboring suboptimal control as a guidance law is discussed 
in terms of the lunar launch problem which is defined in the next section. 
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LUNAR LAUNCH PROBLEM 

The lunar launch problem is to insert a payload in circular lunar orbit over a 
flat moon using a rocket with constant thrust acceleration. While this is a free tina 
time problem, it can be converted to a “fixed final time” problem by choosing the 
horizontal component of velocity as the variable of integration. With the variable 
of integration normalized as « = (« - «„)/(«/ - «o), the optimal control problem 
is stated as follows: Find the thrust inclination history 0(u) which minimizes the 
performance index 


subject to the equations of me 

dt_ 

du 

dy_ 

du 

dv 

du 

and the boundary conditions 

u 0 = 


u } - 1, 

In these equations, a = 20.8 ft/sec 2 is the thrust acceleration, g = 5.32 ft /sec 2 1S the 

acceleration of gravity, u } = 5444 ft/sec is the satellite speed, and u 0 - 0 ft/sec. 

For a piecewise linear control involving nine nodes, the nonlinear programming 
code VF02AD gives the following suboptimal control in degrees: 

0! = 26.0 1 0 2 = 23.3 1 0 3 = 20.51 

0 4 = 17.65 0 5 = 14.86 0 6 = H-90 ( 24 ) 

0 7 = 8.98 0 8 = 6.01 0 9 = 3.03 

Two approaches for applying neighboring suboptimal control are discussed: the 
single optimization approach and the multiple optimization approach Here, u 0 - 0 
for the single optimization approach or a node value for the multiple optimization 
approach. In Ref. 4, neighboring suboptimal control results are presented for the cases 
where there is a thrust acceleration or a gravity modeling error. Only t e gravi y 
case is discussed here because it has the largest errors. 


J = tj 

)tion 

(18) 

(Uf - Uo) 

(19) 

a cos 9 

(u f - U 0 )v 

(20) 

a cos 9 

(uj — u o )(asin0 — g) 

(21) 


OCOS0 


0, to = o, yo = 0, u 0 = 0 , 

(22) 

y, = 50,000 ft, vj= 0 ft/ sec. 

(23) 



SINGLE OPTIMIZATION APPROACH 

In this approach, the suboptimal control from node 1 to node 9 is considered to 
be a sequence of reduced-node suboptimal controls. In other words, the suboptimal 
control from node 1 to node 9 is a nine-node suboptimal control. From node 2 to 
node 9, it is an eight- node suboptimal control; from node 3 to node 9, it is a seven- 
node suboptimal control; and so on. At node 8, there are only two nodes available, 
but these are enough to satisfy the boundary conditions (no optimization). Next, the 
9x3 gain matrix, K 0 in Eq. (17), is computed backward to each node and saved. 1 he 
gains associated with the state t are all zero because there is no condition imposed 
on t f . Hence, the gain matrix, reduces to a 9 x 2 matrix, and the states are now 

6x1! = [<$2/o . r t J 

If the state perturbation occurs at node 8, only Sa 8 is of interest for a sample an 
hold system. Hence, only the gains K 0 ( 8, 1) and A' 0 (8,2) are needed. Similarly, if the 
state perturbation occurs at node 7, only K 0 ( 7, 1) and K 0 (7,2) are needed to compute 
Sa 7 , and so on. For a state perturbation between nodes, the gains are obtained by 
linearly interpolating the gains at adjacent nodes. To have gains over the last or 8 
interval, the gains at nodes 7 and 8 are linearly extrapolated. In conclusion, only the 
gains K 0 (i, 1) and K 0 (i, 2) where i = 1, . . . , 9 need to be stored in the flight computer. 

This approach to neighboring extremal control is tested by introducing a T % 
error in the acceleration of gravity. In other words, the true value of g is taken to be 
x5% different than the value being used in the computation of the gains- Gains are 
computed and stored at every node or at every 0.125u for 9 nodes (Table 1). The 
sample points are assumed to occur at every integration step of the simulation. Here, 
64 integration steps are used so that a sample point occurs every 0.015625u. The 
nominal states are obtained by numerical integration of the equations of motion sub- 
ject to the suboptimal control (24). The true states are obtained by integrating the 
equations of motion with the true acceleration of gravity subject to the neigh oring 
suboptimal control. At each sample point, the true states and nominal states are 
differenced and the differences multiplied by the gains to obtain the control pertur- 
bation. The control perturbation is assumed constant over the sample period, but 
it is added to the piecewise-linear nominal control. Hence, the applied control varies 

linearly over the sample period. , . 

The deviations between the true states and the desired values at the final point are 

presented in Table 2 along with the values which would have been obtained had the 
nominal control (24) been applied open loop. On a relative basis, the improvement 
is substantial. However, a statement about the absolute quality of the closed-loop 
results cannot be made without some performance criteria, say for example, that the 
vehicle has only so much AV to meet the desired final conditions precisely. 

Time histories of the deviations are shown in Fig. 1. Throughout the trajectory, 
the deviations are small, but they do not go to zero at the end. There are two 
possible reasons for this: (a) the quality of the suboptimal trajectory as the vehicle 
moves along its path and (b) the size of the last interval over which the gains are 
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Table 1 


9-NODE SINGLE OPTIMIZATION GAINS 


Node 

y Gain 

v Gain 

1 

-0.369E-5 

-0.673E-3 

2 

-0.289E-5 

-0.462E-3 

3 

-0.385E-5 

-0.521 E-3 

4 

-0.573E-5 

-0.640E-3 

5 

-0.940E-5 

-0.831E-3 

6 

-0.179E-4 

-0.118E-2 

7 

-0.461E-4 

-0.201E-2 

8 

-0.267E-3 

-0.581E-2 

9 

-0.488E-3 

-0.961E-2 


obtained by extrapolation. u prir( > 

Both of these concerns can be addressed by racreasing the number ^ nodes He , 

the computations have been repeated for 17 nodes. The final pomt ^^.ons^e 
presented in Table 2 and show considerable improvement relative to 
However, the deviation histories do not change appreciably relative to Fig. • 


Table 2 


DEVIATION FROM DESIRED FINAL CONDITIONS 


Closed Loop 


% Change 
in g 

State 

Open Loop 

9 Node 
Single Opt 

-5.0 

yj 

9891.024 

65.178 


Vf 

72.540 

-2.705 

+5.0 

vs 

-9891.023 

-63.9S9 

Vf 

-72.540 

2.616 


Closed Loop 
17 Node 
Single Opt. 
20.959 
-1.977 
-19.917 
1.832 


Closed Loop 
9 Node 
Mult. Opt. 
48.137 
-1.566 
-47.891 
1.542 
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Figure 1: 9-Node Single Optimization Deviation Histories 


MULTIPLE OPTIMIZATION APPROACH 

In an attempt to improve just the quality of the neighboring suboptimal control a 
9-node suboptimal control to the final constraint manifold is computed from each node 
of the nominal trajectory (Fig. 2), and the gains are computed for each subtrajectory 
by Eq. (17). These gains are presented in Table 3 and are seen to be larger than those 
of the single optimization approach and uniformly increasing toward the final poin . 
The use of these gains in the simulation with a =f 5% mismatch in the acceleration o 
gravity leads to the final results of Table 2. These closed-loop results are somewhat 

better than those of the single optimization results for 9 nodes. 

The time histories of the deviations are shown in Fig. 3. Overall these deviations 
axe smaller than those of single optimization. Again, the fact that the deviations o 
not go to zero can probably be attributed to the extrapolation of the gains at nodes 

7 and 8 over the last interval. 


DISCUSSION AND CONCLUSIONS 

Two approaches for computing the gains for the neighboring suboptimal control 
guidance law have been tested in a simulation of a lunar launch vehicle: the sing e 
optimization approach and the multiple optimization approach. In both approaches, 
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Figure 2: Multiple Optimization Approach 


Table 3 

9-NODE MULTIPLE OPTIMIZATION GAINS 


Node 

y Gain 

v Gain 

i 

-0.369E-5 

-0.673E-3 

2 

-0.494E-5 

-0.780E-3 

3 

-0.688E-5 

-0.921E-3 

4 

-0.101E-4 

-0.112E-2 

5 

-0.161E-4 

-0.141E-2 

6 

-0.290E-4 

-0.190E-2 

7 

-0.661E-4 

-0.288E-2 

8 

-0.267E-3 

-0.581E-2 

9 

-0.468E-3 

-0.874E-2 
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Figure 3: 9-Node Multiple Optimization Deviation Histories 


a suboptimal control and trajectory with evenly spaced nodes is used as a base, and 
the number of gains which must be stored is very small. 

For single optimization, that part of the suboptimal trajectory from a generic 
node to the final constraint manifold is thought of as a reduced-node suboptimal 
trajectory. Hence, the control becomes less optimal (fewer nodes) toward the end of 
the trajectory and eventually runs out of nodes for satisfying the boundary conditions. 
However, the gains generated by this approach produce good results in a guidance 
simulation. The final point results can be improved by increasing the number of 
nodes. 

The multiple optimization approach is to find a full-node suboptimal control from 
each node of the nominal path to the final constraint manifold. Gains generated from 
these subtrajectories are larger than those of the single optimization approach, are 
uniformly increasing toward the final point, and produce better guidance results, that 

is, the deviations are smaller along the path. 

From these results, it is apparent that the single optimization approach can satis- 
factorily meet the final conditions. On the other hand, if the perturbed trajectory is 
to lie close to the nominal trajectory, the quality of the optimization along the path 
must be improved. Multiple optimization does this, but the amount of computation 
is considerably more than that of single optimization. 
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